#include "SQ_ZZ_pEXFactoring.h"


/**
 *	evaluates a polynomial in x at another polynomial, i.e. g(b) mod f
 **/
ZZ_pEX	SQ_eval(const ZZ_pEX& f,const ZZ_pX& g,const ZZ_pEX& b)
{
	ZZ_pEX res,Q;
	for(int i=0;i<=deg(g);i++)
	{
		if( !IsZero(coeff(g,i)) )
		{
			ZZ_pEX temp;
			SQ_divr(SQ_power(b,i),f,Q,temp); 		
//cout << "temp: " << temp << endl;
			res+=coeff(g,i)*temp;
		}
	}	
//cout << "res: " << res << endl;
	return res;
}


/**
 * given a row-reduced echelon form, solves the homogenous system of equations Ax=0 with GF(2) as fixed field
 */
vec_ZZ_p solve_equ_GF2(mat_ZZ_pE A,int* leading_ones_pos)
{
	// solution vector
	vec_ZZ_p b; // b=0
	int m = A.NumRows();
	int n = A.NumCols();
	b.SetLength(n);
	ZZ_p p_ONE;	// p_ONE = 1
	set(p_ONE);
	int mask[n];	// mask to keep track of changes
	memset(mask,-1,sizeof(int)*n);
	
	for(int i=m-1;i>=0;i--)
	{
		if( !IsZero(A[i]) )
		{
			mask[leading_ones_pos[i]]=1;
			for(int j=n-1;j>leading_ones_pos[i];j--)
			{
				if( !IsZero(A[i][j]) )
				{
					if( mask[j]==-1 ) 
					{
						b[j] = p_ONE;
						mask[j]=1;
					}
					b[leading_ones_pos[i]] += p_ONE;		
				}
			}		
		}
//cout << "mask: " << endl;
//for(int i=0;i<n;i++)
//	cout << mask[i] << " ";
//cout << endl;
//cout << "b: " << endl << b << endl;
	}
	for(int i=0;i<n;i++)	// set independent variables to 1
		if(mask[i]==-1)
			b[i] = p_ONE;

	return b;
}


/**
 * given a row-reduced matrix, reduces it to reduce rwo echelon form, A_{mxn}
 */
int* reduced_row_echelon(mat_ZZ_pE& A,int p,int order)
{
	int m = A.NumRows();
	int n = A.NumCols();
	int* leading_ones_index=new int[m];
	for(int i=0;i<m;i++)
	{
		if( !IsZero(A[i]) )	// traverse non-zero rows
		{
			long j=0;
			while( j<=n && IsZero(A[i][j]) ) j++;	// find location of leading coefficient
			//  make leading coefficient a 1
			if( !IsOne(A[i][j]) )
			{
//				ZZ_pX temp_poly;
//				for(int k=0;k<A[i].length();k++)
//					SetCoeff(temp_poly,k,A[i][k]);
//				MakeMonic(temp_poly); 
				A[i]=power(A[i][j],p-2)*A[i]; // x^{p-1}=1
			}
			leading_ones_index[i]=j;	// keep track of location of leading ones in each row for sorting 
			// reduce
			for(int k=0;k<m;k++)
			{
				if(k!=i && !IsZero(A[k][j]) )
					A[k]-= A[k][j]*A[i];		// assuming leading coefficient is 1 !!!!!!!!!!!!!!!!
			}
		}
	}
	if(order)
	{
		for(int i=1;i<m;i++)
			for(int j=0;j<i;j++)
				if(leading_ones_index[i]<leading_ones_index[j])
				{
					int temp_index=leading_ones_index[i];
					vec_ZZ_pE temp_row = A[i];
					leading_ones_index[i]=leading_ones_index[j];
					A[i] = A[j];
					leading_ones_index[j]=temp_index;
					A[j] = temp_row;
				}
	}
	return leading_ones_index;
}

int SQ_Giesbretch_fac(const ZZ_pEX& f,int p,int mu,ZZ_pEX& b1,ZZ_pEX& b2)
{
	int n = (int) deg(f);
	int dim_of_W = n*mu;
	// theta = primitive element
	ZZ_pE theta = random_ZZ_pE();
	while( IsZero(theta) || IsOne(theta) )
		theta = random_ZZ_pE();
	theta = power(theta,(long) p*mu);	// non-efficient way of geerating the primitive element using x^{p*mu}=theta, x!=0,1
//cout << endl << theta << endl;
	// create extension field F with primitive element \theta
	ZZ_pEX W[dim_of_W];
	for(int i=0;i<n;i++)
	{
		for(int j=0;j<mu;j++)
		{
			W[i*mu+j] = ZZ_pEX(i,power(theta,j)); // \theta^{j}x^{i}
//cout << W[i*mu+j] << endl;
		}
	}

	// Form a basis for K-algebra (K is the fixed field of F under automorphism)
	// find tranformation matrix T such that T(u)= 0 mod f for every basis element of F[x,theta]/(f)
	ZZ_pEX T_u,Q;
	mat_ZZ_p T;
	T.SetDims(dim_of_W,dim_of_W);
	for(int i=0;i<dim_of_W;i++)
	{
		SQ_divr(SQ_mulr(f,W[i]),f,Q,T_u);
		for(int j=0;j<n;j++)
		{
			ZZ_pE pE_coeff = coeff(T_u,j);	// inefficient  way of getting coefficients in NTL library
			ZZ_pX p_coeff =  rep(pE_coeff);
			for(int k=0;k<mu;k++)	// set rows of T with transformed basis vectors
				T[i][j*mu+k]= coeff(p_coeff,k);
		}
	}
//cout << "T: " << endl << T << endl;
	// find kernel of T
	mat_ZZ_p T_kernel;
	kernel(T_kernel,T);	
//cout << "kernel: " << endl << T_kernel << endl;


	// Find a non-zero zero-divisor in K-algebra
	int dim_of_Kalgebra = T_kernel.NumRows();

	// generate random a_{1} and a_{2} in K-algebra
/************************************************
*
*	Random Coefficient Generation
*************************************************/
//cout << "Enter rand vector of size " << dim_of_Kalgebra << ": ";
vec_ZZ_p rand_poly_coeff;
rand_poly_coeff.SetLength(dim_of_Kalgebra);
vec_ZZ_p x;
ZZ_pX random_poly;
random(random_poly,dim_of_Kalgebra);
for(int i=0;i<dim_of_Kalgebra;i++)			// more efficient way to assign poly coefficients to a vector ???!??!!!!
	rand_poly_coeff[i] = coeff(random_poly,i);	// 	?!?!? make sure constant coeff is non-zero ?!?!???????????!!!!!!!!!
//cin >> rand_poly_coeff;	
	mul(x,rand_poly_coeff,T_kernel);
cout << "x: " << x << endl;
	ZZ_pEX a[3];
	for(int i=0;i < T_kernel.NumCols(); i++)
		a[0]+= x[i]*W[i];
cout << "a1: " << a[0] << endl;

random(random_poly,dim_of_Kalgebra);
for(int i=0;i<dim_of_Kalgebra;i++)			// more efficient way to assign poly coefficients to a vector ???!??!!!!
	rand_poly_coeff[i] = coeff(random_poly,i);	// 	?!?!? make sure constant coeff is non-zero ?!?!???????????!!!!!!!!!
//cin >> rand_poly_coeff;	
	mul(x,rand_poly_coeff,T_kernel);
cout << "x: " << x << endl;

	for(int i=0;i < T_kernel.NumCols(); i++)
		a[1]+= x[i]*W[i];
cout << "a2: " << a[1] << endl;

a[2]=SQ_mulr(a[1],a[0])-SQ_mulr(a[0],a[1]);
cout << "a2a1-a1a2: " << a[2] << endl;
/*************************************************/


for(int i=0;i<3;i++)
{
   ZZ_pEX b;
	// find first nu powers of b mod f where  b \in {a_{1},a_{2},a_{1}a_{2}-a_{2}a_{1}}/{0}, these powers of linearly dependent for coefficients \alpha_{0},\alpha_{1},...,\alpha_{nu}
   if(!IsZero(a[i]))
   {
	b = a[i];
	ZZ_pEX pows_of_b[dim_of_Kalgebra+1];
	pows_of_b[0] = SQ_power(b,0);
	// find b^{1},b^{1},...,b^{nu}
	for(int i=1;i<=dim_of_Kalgebra;i++)
	{
		SQ_divr(SQ_mulr(b,pows_of_b[i-1]),f,Q,pows_of_b[i]); 
cout << "b^{" << i << "}: " << pows_of_b[i] << endl;	
	}
      	// solve lin dep equation for \alpha_{0},\alpha_{1},...,\alpha_{nu}
	vec_ZZ_p lin_combo;
	int found_min_poly=0;
	for(int i=1;i<=dim_of_Kalgebra;i++)
	{
		mat_ZZ_pE A;
		A.SetDims(i+1,n);
		vec_vec_ZZ_pE vec_equ_system;
		for(int j=0;j<=i;j++)		// optimize: no need to make the matrix every time !!!
		{
			vec_ZZ_pE vec_equ_system_temp= VectorCopy(pows_of_b[j],n);	// use first i powers of b
			append(vec_equ_system,vec_equ_system_temp);
		}
		MakeMatrix(A,vec_equ_system);  // set up the system of equation: Ax=0
		A = transpose(A);
cout << "A: " << endl << A << endl;
		long rank = gauss(A);
		// make each row as a polynomial monic
		

		// find solutions
cout << "gauss A: " << endl << A << endl;
		int* leading_ones = reduced_row_echelon(A,4,1);
cout << "reduced A: " << endl << A << endl;
		lin_combo = solve_equ_GF2(A,leading_ones);
cout << "solution: " << endl << lin_combo << endl;
		if( !IsZero(lin_combo) )
		{
			found_min_poly = 1;
			break;
		}
	}
	// factor minimal poly
	if(found_min_poly)
	{
cout << "Found minimal polynomial ..." << endl;
		ZZ_pX min_poly;
		for(int i=0;i<lin_combo.length();i++)
			SetCoeff(min_poly,i,lin_combo[i]);
		vec_pair_ZZ_pX_long factors;
		CanZass(factors,min_poly);
if(min_poly == factors[0].a) // if f is irreducible
{
	if( deg(min_poly) == dim_of_Kalgebra )
	{
		cout << "!!! K-algebra is a field !!!" << endl;
		return 0;
	}
}
else	// if f is reducible
{
	int num_of_factors = factors.length();
for(int i=0;i<num_of_factors;i++)
{
	cout << "factor " << i << ": " << factors[i] << endl; 
	ZZ_pEX fac = SQ_eval(f,factors[i].a,b);
	cout << "f factor: " << fac << endl;
	MakeMonic(fac);
	cout << "monic f factor: " << fac << endl;
	
ZZ_pEX Q,R;
	SQ_divr(f,fac,Q,R);
cout << f << " = (" << Q << ")(" << fac << ") + (" << R << ")" << endl;
}
	// f = gh where b1=g, b2=h
	ZZ_p p_ONE;
	set(p_ONE); // p_ONE = 1 
	ZZ_pX g=ZZ_pX(0,p_ONE);
	ZZ_pX h=ZZ_pX(0,p_ONE);
	for(int i=0;i<num_of_factors/2;i++)
		g*=factors[i].a; //power(factors[i].a,factors[i].b);		// !!???!! not taking into acount the exponents
cout << "g: " << g << endl;
	for(int i=num_of_factors/2;i<num_of_factors;i++)
		h*=factors[i].a; //power(factors[i].a,factors[i].b);
cout << "h: " << h << endl;
//check
if( min_poly != (g*h) )	cout << "!!! ERROR: min_poly != g*h" << endl;

	b1 = SQ_eval(f,g,b);
	b2 = SQ_eval(f,h,b);
cout << "b1: " << b1 << endl;
cout << "b2: " << b1 << endl;
	return num_of_factors;
}
	}
	else
		cout << "No minimal polynomial was found!!!" << endl;

   }
}
	cout << "!!! Failure !!!" << endl;
	return -1;
}
